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Using the Green's function method, we study the effect 
of an impurity potential on the electronic structure of the 
honeycomb lattice in the one-band tight-binding model 
that contains both the nearest neighbor (t) and the sec- 
ond neighbor {t') interactions. The model is relevant to 
the case of the substitutional vacancy in graphene. If the 
second neighbor interaction is large enough (t' > t/3), 
then the linear Dirac bands no longer occur at the Fermi 
energy and the electronic structure is therefore funda- 
mentally changed. With only the nearest neighbor inter- 
actions present, there is particle-hole symmetry, as a re- 
sult of which the vacancy induces a "zero-mode" state 



1 Introduction There is considerable current interest 
on the behavior of graphene, which is a two-dimensional 
honeycomb lattice of carbon atoms. The honeycomb lat- 
tice results in the unusual linearly-dispersive band struc- 
ture, which leads to a host of interesting properties directly 
originating from the linearity of the band structure such as 
Klein tunnelling, Zitterbewegung, novel Quantum Hall ef- 
fect, and tabletop quantum electrodynamics |!T!l2l. In this 
paper, we study the impurity states, which also shows un- 
usual behaviors originating from the linearity of the band 
structure and the particle-hole symmetry present in the bi- 
partite honeycomb lattice. Since the physical properties of 
a material are often controlled by the presence of impuri- 
ties, the study of the single vacancy forms the basic founda- 
tion for the understanding of the behavior of more complex 
defects. 



at the band center with its wave function entirely on the 
majority sublattice, i. e., on the sublattice not containing 
the vacancy. With the introduction of the second neigh- 
bor interaction, the zero-mode state broadens into a res- 
onance peak and its wave function spreads into both sub- 
lattices, as may be argued from the Lippmann-Schwinger 
equation. The zero-mode state disappears entirely for the 
triangular lattice and if t' is large for the honeycomb lat- 
tice as well. In case of graphene, t' is relatively small, so 
that a well-defined zero-mode state occurs in the vicinity 
of the band center 
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There have been several theoretical studies of the va- 
cancy center in the honeycomb lattice using the tight- 
binding model |l3]|4l|5l. But, they are either restricted to 
the nearest neighbor (NN) interactions only or they use 
finite-size supercell calculations. In the latter case, it be- 
comes difficult to separate the effects of the finite size from 
the effects of the presence of the impurity. One interest- 
ing property, however, has been clearly established in the 
form of a theorem, viz., the presence of the zero-mode 
state, which states that if only NN interactions are present 
in a bipartite lattice, then the substitutional vacancy pro- 
duces a "zero-mode" state with energy equal to the on-site 
energy, with its wave function living entirely on the sublat- 
tice opposite to that of the vacancy |[3]|5|. In addition, the 
zero-mode state is quasilocalized with the wave function 
falling off as 1/r with distance, in the approximation of the 
linear band structure. However, the effect of the absence 
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Figure 1 The honeycomb lattice with two different sub- 
lattices, A and B, shown as full and open circles and the 
corresponding Brillouin zone with the Dirac points K and 
K' . Impurity is placed at the site Aq and two of its neigh- 
boring sites are labelled as Ai and Bi. 



of the particle-hole symmetry, which is the case if higher- 
neighbor interactions are present, on the zero-mode state 
has not been adequately studied. In this paper, we study in 
detail the effect of the second neighbor {2NN) interaction 
on the zero-mode state. If the 2NN interaction is strong, 
then the honeycomb lattice shows the behavior of the tri- 
angular lattice, where the zero-mode state disappears. The 
presence of the 2NN interaction breaks the particle-hole 
symmetry in a bipartite lattice and the inclusion of the 
further neighbor interactions do not change the physics in 
any essential ways. 

The method we will use in this paper is the Green's 
function (GF) approach, which has been extensively used 
for defects in solids and has been pioneered by Sankey and 
coworkers ||6] for the study of defects in semiconductors. 

2 Tight-binding electronic structure Honeycomb 
lattice is a two-dimensional structure formed by two inter- 
penetrating triangular lattices, A and B, with two atoms in 
the unit cell indicated in Fig. ([1]). The tight-binding Hamil- 
tonian is 



NN 



2NN 



(1) 



keeping only the NN and the 2NN terms, where c^^ , Ci^ 
are the creation and annihilation operators for the elec- 
trons with ia being the cell and sublattice index. From 
density-functional calculations, the hopping parameters for 
graphene are t = 2.56 eV and t' = 0.160 eV |3. 

The band structure is simple to obtain. Using the 
Bloch function basis \ka) — N^^^^ gik.ria jj^,^^ where 
Via = Ri + Ta are the atom positions and TV is the number 
of unit cells in the crystal, the Hamiltonian becomes 



Ho — 



f'ik) f{k) 

r{k) f{k) 



(2) 



with f{k) = -t Y^^j^^ e^'^ '^J = -t[2exp{ikya/2)cos 
(V3fc^a/2)+cxp(-ifcya)] and/'(fc) = t' e^'^ '^i 
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Figure 2 The band structure of the honeycomb lattice for 
different values of the hopping parameters. For t' > t/3, 
the linearly-dispersive Dirac band structure occurring at 
Ep for the half-filled system is lost. With t' /t oo the 
honeycomb lattice becomes effectively two decoupled tri- 
angular sublattices, whose band structure is shown by the 
dash-dotted line. The Fermi energy is indicated in the fig- 
ure for the half-filled system (one electron per site). 

t' [1 cos{\/ikxa) + 4 cos(>/3fc2;a/2) cos(3fcj,a/2)], where 
di's are the NN positions for f{k) and the 2A^A^ posi- 
tions for f'{k). These two functions happen to be related: 
f'{k) = t'|/(fc)p/i^—3t', so that the eigenstates are given 
by 

3t' 



e±ik) 



V2 



±i/(fc)h 

1 / ±e' 



1 



(3) 



where ± denotes the conduction and the valence bands and 
the phase factor 6'^*= = f{k)/\f{k)\ justifies the interpre- 
tation of the wave function in terms of pseudo-spins. 

The band structure is shown in Fig. (|2}. The linear- 
dispersion near the Dirac points K or K', viz., e{q) = 
ztvpq — 3t', apart from the unimportant energy shift, is 
clearly seen in the band structure and the Fermi velocity 
vp = 3ta/2 does not depend on t'. However, if the 2NN 
interaction is strong, the bands cross the zero of energy at 
points other than the Dirac points thus moving the Fermi 
energy Ep away from the linearly-dispersive Dirac points, 
so that the interesting physics originating from the Dirac 
points is lost. This happens when t' > i/3, so that in 
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that case the band structure of the honeycomb lattice be- 
comes fundamentally different, since there is no linearity 
of the band structure at the Fermi energy. The parameters 
for graphene are far from this condition, so that the Dirac 
point physics continues to remain valid even under exter- 
nal perturbations such as strain or electric fields. For very 
large values of the 2NN interaction t << t' , one recov- 
ers the band structure of the triangular lattice as the two 
sublattices become effectively decoupled. 

3 Green's functions and the impurity states The 

impurity is modeled by adding an on-site perturbation V to 
the Hamiltonian 

n^Ho + V, (4) 

where V = C/qcJ^coa, with the impurity located in the 
central cell on the A sublattice and Uo being the strength of 
the impurity potential (vacancy corresponds to Uq oo). 

The key quantity to compute is the full GF G, which 
is related to the unperturbed GF G" through the Dyson's 
equation G = G° + G^VG. With the localized impurity 
potential, the Dyson's equation becomes 



(5) 



{ia\G^ {E)\j P) . This equation may 
be inverted as usual to yield the full GF 



where Gt^^p{E) 



G„ 



1 - C/oGO 



(6) 



To calculate the real-space unperturbed GF, we first 
express it in terms of the momentum GF G'^p{k,E) = 
{ka\G'^{E)\k/3). Using the expression 



G'>iE) = {E + irj-no)-' = J2 



ks 



E 



Ho ' 



(7) 



where s = ± is the band index, and taking the real-space 
matrix elements, we find 

(*a|G°(i?)|j/3) 
1 



n 



BZ 



-^^^^Glp{k,E), 



(8) 



where 
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Figure 3 Impurity state from the Dyson's equation. 
Dashed and solid lines indicate respectively the real and 
imaginary parts - the latter multiplied by — tt^^ to yield 
the DOS po{E) - of the unperturbed GF for three differ- 
ent cases: (a) t ^ l,t' = (top), (b) t ^ l,t' = 0.25 
(middle), and (c) t ~ 0,t' = 1 (bottom). The last case 
corresponds to two decoupled triangular lattices. In the top 
two figures, black dots show the energy of the impurity 
state with ?7o/N = 5 (top) and L/q = 00 (middle). For 
Uq — 00, corresponding to the vacancy, the impurity state 
occurs exactly at Ep in the top figure due to particle-hole 
symmetry, while in the middle figure it is displaced from 
Ep due to lack of symmetry, and in the bottom figure the 
corresponding impurity state has disappeared. 



GlJk,E)=t,k,E 



E 



zr]-f{k) 
f*{k) E 



f{k) 

-tV-fik) 



(9) 



with ^k,E = [{E + ir]- f'{k)f - |/(fc)|2]-i. Thus the 
calculation of G^^ boils down to the numerical integra- 
tion in Eq. (|8]i and the full GF with the impurity is obtained 
from Eq. (|6]l. The faster Horiguchi method may be used for 
the evaluation of the unperturbed GF for the case with NN 
hopping only ||8]|9l. 



All quantities of interest may be expressed in terms of 
the full GF, e. g., the local density of states (LDOS) at a 
specific site is expressed as pict(i?) ~ ~TT^^liaGia^ia{E). 
From Eq. (|6]l, the LDOS at the impurity site has an espe- 
cially simple form 



Poa{E) 



PoiE) 



(10) 



{l-UoFoiE))^ + {7TUoPom^' 
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where po{E) = — tt ^Im Gq^q^(£') is the unperturbed 
DOS, which is the same for every site, and Fo{E) = 
Re Gq^ Qy^{E). According to Eqs. (|6]l and ( fTOl i. the con- 
tribution of the impurity has a sharp peak at the resonance 
energy Eo satisfying the resonance condition 



1 - UoFoiEo) = 0. 



(11) 



The sublattice-projected densities of states pa{E) or 
Pb{E) may be obtained by taking the trace of the GF: 



'(E) — GiaSajE) 



(12) 




and this can be expressed in terms of the central-site GF 
in real and momentum space, Gq^ oa(-^) ^AAi^^ 
respectively, as shown in our earlier work ||5]- The expres- 
sion for the total DOS is obtained from the sum p{E) = 
Pa{E)+pb{E) to yield 

^ , 1 -Uo[UopoiE)F^{E) + p'„{E){l - UpFojE)) 

(13) 

where the prime denotes the derivative. 

The graphical solution of Eq. ( fTTT i for several cases are 
shown in Fig. (O. Solution Eq determines the resonance 
(inside the continuum band) and bound states (outside the 
continuum band) in the presence of the impurity. Accord- 
ing to Eq. ([Tol l, the solutions at the van Hove singularities 
will produce a negligible change in the LDOS, since ei- 
ther the real or the imaginary part of the unperturbed GF 
diverges at these points. 

Several points may be made from Fig. Q, where we 
have shown the evolution of the GF as we go from the hon- 
eycomb to the triangular lattice by changing the strength of 
t': 

(i) The top figure, which shows the solution of the 
Dyson's equation with only the NN interaction, shows that 
as the impurity potential Uq is increased to infinity, the so- 
lution Eq indicated by the black dot approaches the zero 
of energy resulting in the so-called "zero-mode" state. The 
particle-hole symmetry demands that the real part of the 
GF Fo{E) is antisymmetric about the zero of energy, so 
that the resonance state occurs at E^ = just from symme- 
try. Because the unperturbed DOS po (E) is zero at the res- 
onance, it results in a sharp resonance state of zero width 
as seen from Fig. |4] (top panel). 

(ii) With the second 2NN interaction present, particle- 
hole symmetry is no longer present, resulting in a DOS 
that is no longer symmetric about Ep, so that the real part 
of the GF is not required to be antisymmetric. This can be 
seen from the relation Re G(r,£;) = Tr^ip J'^^dE'{E'- 
£')^^Im G(r, E'). This is seen from the middle figure in 
Fig.[3j where the solution Eq occurs no longer at the band 
center, but is shifted from it, the sign of the shift depending 
on the relative signs of t and t'. The resonance peak now 
acquires a finite width being in the band continuum as seen 
from the middle part of Fig. |4] 
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Figure 4 Sublattice densities of states (black solid and red 
dashed lines) together with the unperturbed DOS (black 
dashed line). With the presence of the INN interaction, 
the localized zero-mode state (top) turns into a finite-width 
resonance state (middle), which disappears if t' is large 
(bottom). The impurity potential was taken to be Uajt = 
100 in the figure. 



(iii) As discussed in the previous section, if i' > t/3, 
the band structure is drastically changed with the Fermi 
energy no longer occurring at the Dirac points, even though 
the bands may be linear there. For somewhat larger value 
of t', the zero-mode solution at the band center disappears 
altogether. This may be seen from the bottom parts of Figs. 
(|3]l and (IDi, which show the GF and the DOS, respectively, 
where t' is substantially large. 

As shown by Pereira et al. ||3], the zero-mode state re- 
sides on the B sublattice alone, if only the NN interac- 
tion is present. This is clearly seen from the top part of 
Fig- SI and will be further discussed later when we discuss 
the Lippmann-Schwinger wave function. The emergence 
of the zero-mode state as the impurity potential is gradu- 
ally increased to infinity is shown in Fig. (|5]l. 
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Figure 5 Evolution of the total DOS piE) with the 
strength of the impurity potential; Uo/t ~ (black dashed 
line), 5 (black solid line), and 100 (red dashed line), as ob- 
tained from Eq. (fT3T l {N = 20 was used in the figures). 



Near resonance, the total DOS may be written as a Lon- 
rentzian 



PiE) = 2po{E) 



1 



r 



ttN {E-Eo)^ + r^' 



(14) 



by expanding the real part of the GF: Fq{E) ~ Uq^ + 
Fq{Eo){E — Eo) + valid near the resonance energy. The 
width of the resonance peak is 7^ = —itpo{Eq)/Fq{Eo), 
which is non-zero if the 2NN interaction is present. The 
energy and width of the resonance peak is shown in Fig. 
|7] as a function of t'. The resonance peak disapperas for 
t' >^ 0.4 t or so. 

For the triangular lattice, as seen from the zeros of the 
real part of the on-site GF in Fig. (|3), there is a resonance 
state at the van Hove singularity at Eq = —2t' and there is 
a bound state above the bands at the energy i?o oo as 
Uq ^ oo (the latter state is not shown in our figures). The 
unperturbed and perturbed DOS are shown in Fig. (0. 

4 Wave function of the impurity state As already 
mentioned, if only the NN interaction is present in a bipar- 
tite lattice, a single vacancy in one sublattice (the minority 
sublattice now since there is one atom less) (a) introduces 
a state at the zero of energy (the onsite energy in the tight- 
binding model) and (b) this state furthermore lives on the 
majority sublattice alone. The first result follows from the 
particle-hole symmetry, which ensures the real part of the 
GF to be antisymmetric about the band center E = 0, so 




Figure 6 Unperturbed (red dashed line) and perturbed 
(black solid line) DOS for a vacancy in the triangular lat- 
tice. 




Figure 7 Energy of the resonance impurity state as mea- 
sured with respect to the Fermi energy and the width of the 
resonance peak as obtained from the Lorentzian expression 
Eq. ( fT4l i. All quantities are in units of t. 



that the infinite vacancy potential yields a state at that en- 
ergy following the resonance condition obtained from the 
Dyson's equation. The second result naturally follows by 
examining the perturbed wave function in the presence of 
the impurity potential following the Lippmann-Schwinger 
equation. From this equation, we also estimate the spread 
of the wave function to the minority sublattice if the 2NN 
interactions are present. 

The Lippmann-Schwinger equation relates the per- 
turbed wave function to the unperturbed wave function via 
the GF I If-) = I !?■")+ G^F !!?■). Inverting it, we get 



+ (15) 



and, using the resonance condition Eq. ( fTTT i. obtain the re- 
sult for the resonance state 



ImGO a(Eo)' 



(16) 
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Figure 8 Calculated GF with the 2NN interaction 
present. Real and imaginary parts are shown by dashed and 
solid lines, respectively. Energy scale has been shifted by 
3i' so that the zero coincides with the Fermi energy, which 
continues to occur at the Dirac points for t' < t/3. Results 
for Gq^ were shown in the middle part of Fig. |3] except 
that in that figure the zero of the energy was not shifted. 



To examine the perturbed wave function, we need to ex- 
amine the low-energy behavior of the GF since Eq -> 0. 



Resonance -state wave function for t' = 0: This case 
was discussed in detail in our earlier work ||5]; here we 
briefly summarize the results. For the tight-binding Hamil- 
tonian with only NN interactions, the low-energy GFs are 
given by 

\E\r _ 
\E]ro 



\E\) 



2vf 

r Avp 



i = 0, 



E\E\ 



(17) 



where only the lowest order terms in energy have been 
kept in both the real and the imaginary parts. In Eq. (flTt . 
ro ~ 0.6a, Ac is the unit cell area, and a and /? are numbers 
of the order of one. Plugging this into Eq. (fTSI l and keep- 
ing the dominant terms as Eq — > 0, the perturbed wave 
function becomes 





(18) 



where q is the real part of G^g in Eq. ( [TtI i. Eq. ( fTSb im- 
plies that the vacancy-induced "zero-mode" state resides 
only on the B sublattice and its amplitude is proportional 
to tiie real part of the corresponding GF. Eq. (ITTt also 
contains the striking feature that the wave function decays 
slowly, only as 1 /r with distance, as seen from the real part 

oftheGEGO^^oA- 

Resonance-state wave function for t' ^ 0: We discuss 
the case where t' < t/3, so that the Fermi energy is still 
at the linearly dispersive Dirac point. For higher values 
of t', linearity is destroyed and the resonant-state solution 
near Ep quickly disappears. When the 2NN interaction is 
present, the particle-hole symmetry is lost and the real part 
of the on-site GF is no longer zero at the Fermi energy. The 
imaginary part of the on-site GF is still zero since the band 
structure remains linear with a vanishing DOS at the Fermi 
energy which occurs at —3t'. 

The low energy expressions for the behavior of the GF 
are needed to examine the wave function of the resonance 
state. Unfortunately, it is not easy to derive an analytical 
expression for the GF similar to Eq. (ITTt for the full tight- 
binging bands with the 2NN interactions present. Instead, 
we have studied numerically the behavior of the GFs (some 
typical GFs are shown in Fig. [8]l and extracted the low- 
energy behavior, which goes as 

boE + ico\E\, 



G° 



GO 



iA,QA{E) 
B,OAiE) 



ao 

a^E + ibi\E\. 
di + ie,;i?, 



(19) 

where only the lowest order terms have been kept in both 
the real and the imaginary parts and energy is measured 
from the Dirac point energy —it' , which is the same as Ep 
also. 

From the low-energy behavior of the GF, Eq. ([19), 
we can infer that the resonance state wave function 
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spreads into both sublattices now. To see this, we in- 
sert Eq. ( fT9] l into the Lippmann-Schwinger Eq. (fT6t 
and using the unperturbed wave function Eq. (|3]l, we 
get the result for the perturbed wave function: 'FiA = 



0.5 



(2iV) 



ik-riA 



+ icQ ^{tti + ibi)] and 



(27v)-i/2[eifc-'-.B + ie''^>'CQ^{di/Eo + iet)]. For small t', 
the resonant energy Eq cx t' (see Fig.|7]), so that the dom- 
inant part of the wave function is the second term of 'FiB . 
This leads to the result "FIaI'^iB ~ cld~'^El - {f /tf, 
where we have used the rough behavior cq ^ \/t^ and 
di ^ 1/t. Thus the wave function spreads to the minor- 
ity A sublattice in the presence of the 2NN interaction 
roughly as 





(20) 



This is consistent with Eq. (fTsT l that the wave function re- 
sides on the majority B sublattice alone, when the 2NN 
interaction is absent. 

In Fig. (|9]i, we show the site-specific local DOS for 
sites near the vacancy with the sites labelled in Fig. ([T]i. The 
figure indicates that the wave function for the zero-mode 
state has contributions only from the majority B sublattice 
if = and there is contribution from both the sublattices 
iff 7^ 0. 

5 Summary and Discussions In summary, we stud- 
ied the formation of the impurity states in the honeycomb 
lattice, where both the NN and the 2A^A^ interactions are 
present. In the limit of large 2NN interaction, the honey- 
comb lattice effectively goes over to two decoupled trian- 
gular lattice. If the 2NN interaction is absent, the impu- 
rity in the vacancy limit (impurity potential oo) leads 
to the well-known zero-mode resonance state, introduc- 
ing a (5-function peak in the DOS. As the 2A^A^ interac- 
tion is increased, the resonance state turns into a broadened 
peak, but quickly disappearing altogether for t' >~ {)At. 
For t' = 0, the resonance state is completely confined 
to the majority B sublattice, i. e., opposite to the sublat- 
tice in which the vacancy is present, while with t' ^ 0, 
the resonance-state wave function spreads into both sub- 
lattices. 

For the case of graphene, the 2NN interaction is small 
with t' ^ 0.06t, so that the zero-mode state occurs as a 
sharp peak with energy just above the band center energy 
E = Q. An estimate using the actual hopping parame- 
ters given earlier and results of Fig. (|7) indicates the peak 
position to be _Eo ~ 0.11 eV and a similar peak width 
r « 0.08 eV. In graphene, there are also the sp^a dan- 
gling bond electrons on the three carbon atoms adjacent to 
the vacancy. Three electrons occupy these states and the 
sole remaining electron occupies the zero-mode state dis- 
cussed here. The strength of the Hund's coupling between 
the a and the tt electrons is about 0.4 eV, which would 
pull down the energy of the zero-mode state by 0.2 eV 
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Figure 9 Local DOS at individual sites labelled in Fig. ([T]) 
without (top) and with (bottom) the 2NN hopping. Unper- 
turbed DOS po{E) is shown by the black dashed lines and 
Uq was taken to be lOOi. Top figure indicates that there is 
no contribution to the zero-mode state from the A sublat- 
tice, while in the bottom figure, the Ai site does have a 
small contribution (shown by the arrow). Clearly, the zero- 
mode state which lives on the B sublattice, if only the NN 
hopping is present, has a strong component from the Bi 
site in both cases shown. The contribution of the Aq site is 
of course zero in the limit Uq = oo. 



or so. This means that the zero-mode state becomes spin- 
polarized with one spin state occupied and the other empty, 
leading to a magnetic center Density-functional calcula- 
tions show that this is indeed the case and the vacancy pro- 
duces a magnetic center with the vacancy-induced zero- 
mode TT band state half filled. Recently, experimental evi- 
dence using scanning tunneling microscopy has been ob- 
tained for the existence of the zero-mode state [lOj. We 
also note that irrespective of the position of the resonance 
state, which occurs never too far away from Ep, it is un- 
likely for it to be filled by two electrons considering the 
Coulomb interaction, so that the resonance state is filled 
by one electron and remains magnetic. 

Acknowledgements This work was supported by the 
U. S. Department of Energy through Grant No. DOE-FG02- 
00ER45818. 

References 

[1] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. 
Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009) 



Copyright line will be provided by the publisher 



8 



M. Sherafati et al.: Impurity states on honeycomb lattice 



[2] D. S. L. Abergel, V. Apalkov, J. Berashevich, K. Ziegler, and 

T. Chakraborty, Adv. Phys. 59, 261 (2010) 
[3] V. M. Pereira, J. M. B. Lopes dos Santos, and A. H. Castro 

Neto, Ptiys. Rev. B 77, 115109 (2008) 
[4] M. Hjort and S. Stafstrom, Phys. Rev. B 61, 14089 (2000) 
[5] B. R. K. Nanda, M. Sherafati, Z. Popovic, and S. Satpathy, 

Phys. Rev. B (to be submitted) 
[6] O. F. Sankey and J. D. Dow, Phys. Rev. B 27, 7641 (1983); 

O.F. Sankey, J. D. Dow, and K. Hess, Appl. Phys. Lett. 41, 

664 (1982) 

[7] B. R. K. Nanda and S. Satpathy, Phys. Rev. B 80, 165430 
(2009) 

[8] T. Horiguchi, J. Math. Phys. 13, 1411 (1972) 
[9] M. Sherafati and S. Satpathy, Phys. Rev. B (201 1, in press) 
[10] M. M. Ugeda, I. Brihuega, F. Guinea, and J. M. Gomez- 
Rodriguez, Phys. Rev. Lett. 104, 096804 (2010) 



Copyright line will be provided by the publisher 



